### Power Analysis ###
outcome <- ESS10$mig
run <- ESS10$date
power1 <- rdpower(data=cbind(outcome,run), samph=cbind(3,3))
power2 <- rdpower(data=cbind(outcome,run), samph=cbind(7,7))
power3 <- rdpower(data=cbind(outcome,run), samph=cbind(14,14))
power4 <- rdpower(data=cbind(outcome,run), samph=cbind(21,21))




### RDD ###
rd_mig <- RDestimate(mig ~ date | cntry, data = ESS10, cutpoint = 0, bw = 14, se.type = "HC1")
summary(rd_mig)
rd_eco <- RDestimate(eco ~ date | cntry, data = ESS10, cutpoint = 0, bw = 14, se.type = "HC1")
summary(rd_eco)
rd_cul <- RDestimate(cul ~ date | cntry, data = ESS10, cutpoint = 0, bw = 14, se.type = "HC1")
summary(rd_cul)


### RDD with factor outcome ###
### create factor ###
fact_out <- prcomp(~ESS10$mig + ESS10$cul + ESS10$eco)
fact_out
fact_out_1 <- get_pca_ind(fact_out)
fact_out_1$coord
data_fact <- ESS10 %>% 
  filter(mig >-1 & eco > -1 & cul > -1) 
data_fact$factor <- fact_out_1$coord[,1]
summary(data_fact$factor)

# with canned routing
rd_fact <- RDestimate(factor ~ date | cntry, data = data_fact, cutpoint = 0, bw = 14, se.type = "HC1")
summary(rd_fact)

# manually, to get slopes for regression table
smol_fact <- data_fact %>% 
  filter(date > -15 & date < 15)

smol_fact$wt <- kernelwts(smol_fact$date, center = 0, bw = 14)

# no coefs
slope_fact <- lm(factor ~ afg + date + afg*date + cntry, data = subset(smol_fact, wt > 0), weights = wt)
coeftest(slope_fact, vcov. = vcovHC(slope_fact, type = "HC1"))

# coefs
slope_fact_c <- lm(factor ~ afg + date + afg*date + cntry + eisced + polintr + agea + gndr, data = subset(smol_fact, wt > 0), weights = wt)
coeftest(slope_fact_c, vcov. = vcovHC(slope_fact_c, type = "HC1"))





#customized plot function from rdd package
plot.RD <- function(x,gran=400,bins=100,which=1,range,...) {
  frm<-FALSE
  if("frame" %in% names(x$call)) frm<-eval.parent(x$call$frame)
  if(!frm){
    x$call$frame<-TRUE
    x$call$verbose<-FALSE
    x<-eval.parent(x$call)
  }
  d<-as.data.frame(x$frame)
  
  if(length(x$na.action)>0)
    d<-d[-x$na.action,]
  
  if("kernel" %in% names(x$call)) 
    kern<-eval.parent(x$call$kernel)
  else 
    kern<-"triangular"
  
  if("cutpoint" %in% names(x$call)) 
    cut<-eval.parent(x$call$cutpoint)
  else
    cut<-0
  
  bw<-x$bw[1]
  
  if(missing(range)) {
    range<-c(cut-10*bw,cut+10*bw)
    if(range[1]<min(d$X)) range[1]<-min(d$X)
    if(range[2]>max(d$X)) range[2]<-max(d$X)
  }
  
  if(range[1]=="min")
    range[1]<-min(d$X)
  if(range[2]=="max")
    range[2]<-max(d$X)
  range<-as.double(range)
  
  rdplot<-function(d) {
    d.l<-data.frame(X=d$X[d$X<cut],Y=d$Y[d$X<cut])
    lval<-seq(range[1],cut,length.out=(gran%/%2))
    lest<-vector(length=(gran%/%2))
    llwr<-vector(length=(gran%/%2))
    lupr<-vector(length=(gran%/%2))
    for(i in 1:(gran%/%2)) {
      sub<-d.l$X>=(lval[i]-bw) & d.l$X<=(lval[i]+bw)
      w<-kernelwts(X=d.l$X[sub],center=lval[i],bw=bw,kernel=kern)
      ly<-d.l$Y[sub]
      lx<-d.l$X[sub]
      if(length(lx)<=2)
        pred<-rep(NA,3)
      else
        pred<-predict(lm(ly~lx,weights=w),interval="confidence",newdata=data.frame(lx=lval[i]))
      lest[i]<-pred[1]
      llwr[i]<-pred[2]
      lupr[i]<-pred[3]
    }
    
    d.r<-data.frame(X=d$X[d$X>=cut],Y=d$Y[d$X>=cut])
    rval<-seq(cut,range[2],length.out=(gran%/%2))
    rest<-vector(length=(gran%/%2))
    rlwr<-vector(length=(gran%/%2))
    rupr<-vector(length=(gran%/%2))
    for(i in 1:(gran%/%2)) {
      sub<-d.r$X>=(rval[i]-bw) & d.r$X<=(rval[i]+bw)
      w<-kernelwts(X=d.r$X[sub],center=rval[i],bw=bw,kernel=kern)
      ry<-d.r$Y[sub]
      rx<-d.r$X[sub]
      if(length(rx)<=2)
        pred<-rep(NA,3)
      else
        pred<-predict(lm(ry~rx,weights=w),interval="confidence",newdata=data.frame(rx=rval[i]))
      rest[i]<-pred[1]
      rlwr[i]<-pred[2]
      rupr[i]<-pred[3]
    }
    
    #plot to the left
    if(length(unique(d$Y))==2) {
      #DO THIS for when the outcome is dichotomous
      ep<-(max(d$X)-min(d$X))/(2*bins)
      nX<-seq(min(d$X)-ep,max(d$X)+ep,length=bins+1)
      nY<-rep(NA,length(nX))
      for(i in (1:(length(nX)-1))){
        if(sum(!is.na(d$Y[d$X>nX[i] & d$X<=nX[i+1]]))==0)
          next
        nY[i]<-sum(d$Y[d$X>nX[i] & d$X<=nX[i+1]],na.rm=TRUE)/sum(!is.na(d$Y[d$X>nX[i] & d$X<=nX[i+1]]))
      }
      sub<-nX>=range[1] & nX<=range[2]
      subl<-lval>=range[1] & lval<=range[2]
      subr<-rval>=range[1] & rval<=range[2]
      plot(nX,nY,
           type="p",pch=20,cex=.5,col="black",
           xlim=c(range[1],range[2]),
           ylim=c(min(c(llwr[subl],rlwr[subr]),na.rm=T),
                  max(c(lupr[subl],rupr[subr]),na.rm=T)),
           xlab=NA,
           ylab=NA,
           main=NA
      )
    } else {
      subl<-lval>=range[1] & lval<=range[2]
      subr<-rval>=range[1] & rval<=range[2]
      plot(d$X,d$Y,
           type="p",pch=20,cex=0,col="black",
           xlim=c(range[1],range[2]),
           ylim=c(min(c(llwr[subl],rlwr[subr]),na.rm=T),
                  max(c(lupr[subl],rupr[subr]),na.rm=T)),
           xlab=NA,
           ylab=NA,
           main=NA
      )
    } 
    #plot to the left
    lines(lval,lest,
          lty=1,lwd=6,col="#B22222",type="l"
    )
    
    lines(lval,llwr,
          lty=2,lwd=3,col="#B22222",type="l"
    )
    lines(lval,lupr,
          lty=2,lwd=3,col="#B22222",type="l"
    )
    
    #plot to the right
    lines(rval,rest,
          lty=1,lwd=6,col="#228B22",type="l"
    )
    lines(rval,rlwr,
          lty=2,lwd=3,col="#228B22",type="l"
    )
    lines(rval,rupr,
          lty=2,lwd=3,col="#228B22",type="l"
    )
  }
  if(x$type=="sharp" | 1%in%which){
    rdplot(d)
    dev.flush()
  }
  if(x$type=="fuzzy" & 2%in%which){
    d$Y<-d$Z
    rdplot(d)
    dev.flush()
  }
}

# rdd plot
par(mfrow = c(1,3))
plot(rd_mig, gran = 4, range = c(-14,14))
abline(v=0, col="black", lwd = 4, lty = 1)
plot(rd_eco, gran = 4, range = c(-14,14))
abline(v=0, col="black", lwd = 4, lty = 1)
plot(rd_cul, gran = 4, range = c(-14,14))
abline(v=0, col="black", lwd = 4, lty = 1)




# with covariates
rd_mig <- RDestimate(mig ~ date | cntry + eisced + polintr + agea + gndr, data = ESS10, cutpoint = 0, bw = 14)
summary(rd_mig)
rd_eco <- RDestimate(eco ~ date | cntry + eisced + polintr + agea + gndr, data = ESS10, cutpoint = 0, bw = 14)
summary(rd_eco)
rd_cul <- RDestimate(cul ~ date | cntry + eisced + polintr + agea + gndr, data = ESS10, cutpoint = 0, bw = 14)
summary(rd_cul)

# slopes
smol <- ESS10 %>% 
  filter(date > -15 & date < 15)

smol$wt <- kernelwts(smol$date, center = 0, bw = 14)

slopemig <- lm(mig ~ afg + date + afg*date + cntry, data = subset(smol, wt > 0), weights = wt)
slopeeco <- lm(eco ~ afg + date + afg*date + cntry, data = subset(smol, wt > 0), weights = wt)
slopecul <- lm(cul ~ afg + date + afg*date + cntry, data = subset(smol, wt > 0), weights = wt)

coeftest(slopemig, vcov. = vcovHC(slopemig, type = "HC1"))
coeftest(slopeeco, vcov. = vcovHC(slopeeco, type = "HC1"))
coeftest(slopecul, vcov. = vcovHC(slopecul, type = "HC1"))

slopemigc <- lm(mig ~ afg + date + afg*date + cntry + eisced + polintr + agea + gndr, data = subset(smol, wt > 0), weights = wt)
slopeecoc <- lm(eco ~ afg + date + afg*date + cntry + eisced + polintr + agea + gndr, data = subset(smol, wt > 0), weights = wt)
slopeculc <- lm(cul ~ afg + date + afg*date + cntry + eisced + polintr + agea + gndr, data = subset(smol, wt > 0), weights = wt)

coeftest(slopemigc, vcov. = vcovHC(slopemigc, type = "HC1"))
coeftest(slopeecoc, vcov. = vcovHC(slopeecoc, type = "HC1"))
coeftest(slopeculc, vcov. = vcovHC(slopeculc, type = "HC1"))
